npj Systems Biology and Applications

Constraint-based metabolic modeling is a powerful way to study the mechanistic basis of cellular states and disease, but effective use demands substantial computational expertise and careful coordination of multi-step analyses. We developed MechAInistic to lower this barrier enabling researchers to ask complex biological questions in natural language. MechAInistic is a multi-agent system harnessing large language models organized around an Architect-Reviewer pattern that that converts a natural-language question into an executable, model-grounded workflow and produces a structured report. It supports pathway comparison, perturbation analysis, drug-target exploration, and literature interpretation across healthy and disease paired states. We evaluated MechAInistics therapeutic hypothesis generation using two immune-cell use-cases. For rheumatoid arthritis/healthy Naive B models, it identified mitochondrial metabolic rewiring and nominated Devimistat/CPI-613 as an investigational OGDH-centered hypothesis. In CD4+ Th17 multiple sclerosis/healthy models, the workflow identified NADP-dependent isocitrate dehydrogenase as the optimal target and proposed Ivosidenib as an FDA-approved repurposing candidate. GRAPHICAL ABSTRACT O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=83 SRC="FIGDIR/small/723319v1_ufig1.gif" ALT="Figure 1"> View larger version (19K): org.highwire.dtl.DTLVardef@1b5c1d1org.highwire.dtl.DTLVardef@1c798cforg.highwire.dtl.DTLVardef@10161d3org.highwire.dtl.DTLVardef@1bd7dce_HPS_FORMAT_FIGEXP M_FIG C_FIG

Metabolic networks are typically viewed as homeostatic systems that stabilize flux, energy charge, redox balance, and metabolite availability under perturbation. However, it remains unclear whether the same feedback architectures that support metabolic robustness can also generate learning-like, experience-dependent adaptation. Here, we develop a coarse-grained dynamical model of mammalian energy metabolism to test whether prior perturbation can improve future metabolic responses. The model represents core glucose, glutamine, fatty acid, and oxidative phosphorylation pathways as coupled ordinary differential equations with Michaelis-Menten-type fluxes, product-inhibition feedback, adaptive enzyme-capacity regulation, and explicit ATP costs for enzyme adjustment. Rather than aiming to reproduce quantitative fluxes for a specific cell type, the framework is designed to expose how metabolic feedback, regulatory cost, repeated perturbation, and environmental variability interact. We use this model to ask whether adaptive enzyme regulation enables improved recovery after repeated challenges, whether such effects depend on energetic control costs, and whether environmental variability broadens or constrains the set of reachable adaptive states. This approach provides a tractable way to investigate how homeostatic metabolic regulation may give rise to experience-dependent metabolic plasticity.

Over a thousand fungal genomes have been sequenced, yet manually curated genome-scale metabolic models (GEMs) are available for only a limited number of species. Moreover, these models have often been developed independently, leading to inconsistencies in namespaces, compartment definitions, and pathway representations that hinder comparative analysis, the systematic reuse of prior curation efforts, and the integration of consolidated metabolic knowledge. Here, we present the Consolidated Fungal Core Metabolism Model (CFCMM), constructed by integrating thirteen published fungal models spanning Ascomycota, Mucoromycota, and both Crabtree-positive and Crabtree-negative yeasts. We harmonized metabolites and reactions into a non-redundant shared ModelSEED ontological space, standardized compartmentalization, and refined gene-protein-reaction (GPR) rules. Using pathway-level visualization and systematic gap detection, we further improved the integrated network through literature-guided curation to correct stoichiometry, stereospecificity, and pathway architecture. Orthologous protein family reconstruction and functional annotation workflows were used to validate and inform GPR associations, with particular emphasis on ambiguous enzyme superfamilies and membrane-associated components. Using the resulting CFCMM, we built high-quality central carbon core models for each fungus and performed flux balance analysis to quantify ATP-yield variation under aerobic and anaerobic conditions, explicitly evaluating scenarios driven by differences in electron transport chain (ETC) composition. Simulations reproduced the expected fermentative yield of approximately 2 mmol ATP per mmol glucose under anaerobic conditions and separated the thirteen fungi into two bioenergetic groups under aerobic respiration based on Complex I status, with predicted yields of approximately 30 versus 22 mmol ATP per mmol glucose. Forcing flux through the alternative oxidase bypass further reduced ATP yields to approximately 12 and 4 mmol ATP per mmol glucose in Complex I-containing and Complex I-lacking fungi, respectively. Collectively, this work provides a manually curated, ModelSEED-consistent, and extensible fungal core metabolic template, deployed in DOE KBase as a resource for automated reconstruction of central carbon core models from any sequenced fungal genome. In addition, the CFCMM provides modular components for developing GEMs with more accurate energy predictions and enables robust comparative analyses of fungal bioenergetics and core metabolic diversity.

Protein tyrosine phosphatase receptor J (PTPRJ) restrains cell proliferation and migration by dephosphorylating receptor tyrosine kinases (RTKs) including the epidermal growth factor receptor (EGFR). PTPRJ is a purported tumor suppressor, and alterations to its expression and/or function are associated with colorectal, breast, lung, and other cancers. While there is interest in controlling PTPRJ-regulated phenotypes, efforts are limited by the complexity of PTPRJ-mediated signaling. PTPRJ dephosphorylates multiple RTKs, and the degree to which PTPRJ control of signaling and phenotypes depends on local cellular RTK activation profiles is unknown. To probe the context dependence of PTPRJ signaling regulation, we collected signaling measurements across 16 pathway nodes at two time points in a panel of HSC3 carcinoma cells engineered with different PTPRJ expression profiles. Cells were treated with three different RTK ligands, and paired phenotype measurements (viability, wound healing, xCELLigence cell index) were made. Partial least squares regression models were developed to predict relationships between PTPRJ-regulated signaling pathways and cell phenotypes. The model effectively separated contributions to variance arising from the PTPRJ expression background and growth factor context. In testing model predictions, we demonstrated that PTPRJ suppressed MET-induced cell cell proliferation via regulation of a HER3/AKT signaling axis that stabilized PTPRJ expression through an unanticipated feedback mechanism. We also found that PTPRJ regulated HSC3 cell migration via JNK signaling that was preferentially activated by MET. Our results identify new regulatory nodes through which PTPRJ influences cancer cell phenotypes and demonstrates that these processes preferentially occur in the context of distinct RTK activation states.

Quantitative Systems Pharmacology (QSP) models play an inherently interventional role in pharmaceutical research and development, functioning as executable causal systems for designing, evaluating, and replacing clinical trials. However, deploying QSP as an experimental planning engine remains constrained by the difficulty of translating unstructured literature descriptions of clinical or preclinical scenarios into reproducible, simulation-ready model interventions. Motivated by this issue, we propose an agent-based framework that operationalizes QSP models as intervention-ready experimental systems by automatically extracting and executing literature-derived scenarios. The framework combines semantic grounding of model entities with a large language model (LLM)-driven Scenario Extractor and a dual-agent Scenario Mapper. Rather than relying on opaque, single-shot reasoning, our pipeline converts free-text interventions into precise parameter configurations through discrete, verifiable work orders. Moreover, our dynamic Human-in-the-Loop (HITL) strategy empowers modelers to resolve biological ambiguities interactively. Across four diverse kinetic ordinary differential equation (ODE)/QSP models and seven Subject Matter Expert (SME)-curated literature scenarios, our model resolved all selected scenarios into correct executable parameter changes, including multi-dose interventions, unit conversions, no-op scenarios, and ambiguity-triggered HITL cases, demonstrating that structured collaboration between experts and agentic systems can resolve scenarios that standalone raw Systems Biology Markup Language (SBML) reasoning LLM calls handle unreliably.

Enterococci are Gram-positive opportunistic pathogens responsible for a wide range of nosocomial infections. One enterococcocal species, Enterococcus faecium, is steadily increasing in prevalence and has been listed among major multidrug-resistant ESKAPE pathogens. To gain systems-level insights into its metabolism and support discovery of potential therapeutic targets, we constructed iDR479, a comprehensive manually curated genome-scale metabolic model (GEM) to serve as a digital twin for E. faecium TX0016 (strain DO). The reconstruction was curated through extensive homology searches and literature evidence, and further refined and gap-filled through experimental validation. Phenotypic profiling using Biolog microarrays enabled assessment of carbon source utilization, while amino acid leave-out growth assays allowed the evaluation of auxotrophies. The final refined model is 100% accurate in predicting amino acid auxotrophy and 85% accurate in predicting growth on sole carbon sources. Discrepancies between model predictions and experimental phenotypes identified specific knowledge gaps across metabolic pathways, including unresolved carbon source utilization phenotypes, e.g., psicose, sorbitol, and palatinose utilization. Those gaps will guided future experimental characterization. Additionally, gene essentiality analysis was conducted to evaluate the predictive capacity of iDR479 model. Since no experimental gene essentiality data are currently available for E. faecium, model predictions were compared against Tn-seq experimental results from E. faecalis MMH594. Under simulated rich medium conditions, iDR479 achieved 86.7% concordance with the experimental essentiality results of E. faecalis MMH594. iDR479 thus provides a framework for studying E. faecium, offers insights into its metabolic network, and serves as a source for guiding future research and identification of therapeutic targets.

Microbes play a vital role in plant development, health, and resilience, yet relatively little is known about the specific metabolic mechanisms driving interactions in these host-associated communities. Systems biology models enable a computational approach to understanding metabolic interactions, which can be difficult to pinpoint experimentally; however, these methods cannot yet accommodate the large number of species in natural communities. Synthetic communities (SynComs) provide a more tractable alternative to explore targeted interactions. Here, we investigated metabolite exchange in a seven-member maize root-associated SynCom, specifically accounting for plant host context by designing a customized exudate medium. We constructed metabolic models for each bacterial species and curated them with in vitro phenotyping data to reflect experimentally based carbon uptake potential. Flux balance analysis of individual species demonstrated that integrating phenotype data and changing medium type had substantial impacts on predicted growth rates, which in turn shaped potential interspecies interactions. In silico community growth optimization of the seven-member community model showed that the exudate medium supported a more diverse community composition compared to minimal medium, with predictions of community member abundance closely aligned to literature-derived experimental results. Predicted metabolite exchange in the root exudate environment showed Enterobacter ludwigii as a community hub, and cross-feeding of indole suggested a potential effect of bacterial community interactions on the plant host. Our in silico findings indicate the host plays an important role in structuring microbial interactions and cross-feeding at the metabolic level, underscoring the importance of considering environmental context from both theoretical and experimental perspectives. IMPORTANCETrue understanding of a system is marked by the ability to predict its behavior. The complexity of natural host-microbe systems represents a frontier of knowledge that scientists are working to understand, and elucidating principles of interactions within multi-partite microbial communities remains a challenge in microbial ecology. Synthetic communities provide a tractable starting point for investigating interaction mechanisms, and computational approaches complement laboratory experiments by systematically evaluating multiple possibilities for metabolic pathway processing, thereby allowing us to comprehensively study the interconnected metabolic networks of host-associated microbiota. The model we developed for the seven-member maize root-associated bacterial community presents a step toward predicting plant-microbe behavior, providing hypotheses for future experimental testing and serving as a template for expanding model complexity to more members and other systems.

Intratumoral heterogeneity drives variable drug responses in cancer. Single-cell RNA sequencing (scRNA-seq) enables characterization of such heterogeneity and prediction of drug response at single-cell resolution. Accordingly, various computational models have been developed to infer drug response from scRNA-seq data. However, their performance, robustness, and generalizability across different biological contexts remain insufficiently evaluated. To address this gap, we benchmarked representative single-cell drug response prediction models using 26 curated datasets comprising over 760,000 cells across 12 cancer types and 21 therapeutic agents. We constructed balanced and imbalanced scenarios to reflect realistic drug-response label distributions. To address the lack of ground-truth labels in conventional scRNA-seq datasets, we incorporated lineage-tracing data with experimentally validated drug-response annotations, enabling evaluation in a clinically relevant pre-treatment prediction setting. Our results show that prediction performance was markedly higher in cell lines than in tissue samples. Under imbalanced conditions, most methods exhibited sharp performance declines, whereas scDEAL demonstrated the highest robustness. Independent validation using an in-house pancreatic ductal adenocarcinoma dataset further confirmed scDEALs robustness and ability to capture biologically meaningful state transitions. Label-substitution experiment revealed that this robustness was partially driven by the models specific training-label construction. However, benchmarking with lineage-tracing data revealed a fundamental limitation: most models capture drug-induced transcriptional changes but struggled to predict intrinsic resistance before treatment. In summary, our study defines the performance boundaries of current approaches and highlights their limitations in addressing intratumoral heterogeneity, class imbalance, and intrinsic resistance prediction, emphasizing the need for the next-generation single-cell drug response models with stronger clinical relevance.

Hippocampal adult neurogenesis (HANG) is a highly regulated process where neural stem cells progress through distinct stages--from Type 1 radial glia-like cells to mature neurons--via a complex series of proliferative and differentiative divisions. While recent in vivo imaging has provided valuable insights to cellular processes, the exact relationship between individual cell-fate decisions and long-term population stability remains difficult to quantify empirically. In this study, we utilized an agent-based (AB) model to simulate the stochastic dynamics of the hippocampal neurogenic niche. Our results demonstrate that while individual progenitor lineages exhibit high variability and probabilistic division symmetries (proliferative symmetric, asymmetric, and differentiative symmetric), the system achieves deterministic stability as the initial progenitor density increases. We found that the T1 progenitor pool follows a negative exponential decay profile, with its longevity primarily dictated by the differentiation rate (d,0). Critically, the terminal output of immature neurons (CIN,t) was non-linearly coupled to the proliferative capacity of transit-amplifying cells (pp,0); even marginal increases in symmetric proliferative divisions resulted in an exponential expansion of the neuronal pool. These findings suggest that the homeostatic maintenance of the hippocampal niche is governed by a kinetic tuning of division probabilities, providing a theoretical bridge between single-cell stochasticity and robust tissue-level output.

Cells communicate via extracellular ligands, such as hormones, which bind to plasma membrane receptors and trigger intracellular signaling cascades. G Protein-Coupled Receptors (GPCRs) exemplify this mechanism by initiating signaling both at the cell surface and, from intracellular compartments such as endosomes. The kinetics and spatial localization of these signals are critical determinants of cellular responses, yet receptor trafficking-including internalization, endosomal sorting, and recycling-remains a pivotal but often overlooked component of theoretical GPCR models. In this study, we present a mathematical framework that integrates receptor trafficking and signaling compartmentalization into generic GPCR dynamic models. Using a compartmentalized approach based on systems of ordinary differential equations (Chemical Reaction Networks), we analyze how receptor internalization and recycling modulate ligand-induced responses. Our results show that the balance between plasma membrane and endosomal signaling can significantly enhance or diminish ligand efficacy. Calibrated with high-throughput kinetic data, our model offers a refined tool for ligand pharmacological characterization and advances the understanding of GPCR signaling spatial organization.

Joint modelling of longitudinal data using non-linear mixed effects models and time-to-event outcomes provides a suitable framework to account for informative censoring when estimating biomarker dynamics and quantifying event risk using covariates and longitudinal trajectories. Their usefulness in clinical research depends on data collection design, particularly to precisely estimate the association (link) parameter between longitudinal and survival processes. However, optimal design strategies have so far been addressed separately for longitudinal and survival endpoints and remain unexplored for joint models. We propose two Fisher Information Matrix (FIM) computation methods for joint models, relying on Monte-Carlo integration over observations combined with either Markov Chains Monte-Carlo or Adaptive Gaussian Quadrature to integrate random effects. Their accuracy is assessed against clinical trial simulations in an oncological example based on the HORIZON III study with a tumour-growth-survival model including discrete and continuous covariates. We apply these methods to quantify the impact of follow-up duration, sampling richness, sample size, and covariate distribution on parameter uncertainty and test power. In our example, longitudinal-parameter uncertainty is barely affected by follow-up duration or sampling richness, whereas survival-parameter uncertainty decreases substantially from 1-year to 2-year follow-up. The number of subjects needed (NSN) to achieve <15\% uncertainty on the link parameter is comparable for a 2-year rich design and a 3-year sparse design. Optimal covariate distributions are stable across designs and systematically improve test power, outperforming longer and richer but non-optimised designs. These FIM-based methods accurately predict uncertainty and test powers, enabling design evaluation and NSN computation for joint-model-based clinical studies.

Glioblastoma is an aggressive and highly heterogeneous brain tumor with poor prognosis despite multimodal treatment strategies. Understanding the molecular diversity of the disease is essential for improving tumor stratification and identifying potential therapeutic targets. In this study, we investigate whether network-based analysis can reveal biologically meaningful subgroups of glioblastoma tumors. Using RNA sequencing and mutation data from the TCGA-GBM cohort, we constructed patient-specific protein-protein interaction networks based on genes that are differentially expressed or harbor somatic mutations. These networks capture the molecular alterations associated with individual tumors within the context of the human interactome. We then derived similarities between tumors using a binary representation of network nodes and the Jaccard similarity metric, enabling the construction of a patient similarity graph. Community detection algorithms (Louvain and Leiden) were applied to this graph to identify clusters of tumors with similar molecular network profiles. Our analysis revealed six tumor communities characterized by distinct gene compositions and enriched biological processes. For each community, we identified candidate biomarkers and network hubs that may represent potential therapeutic targets. Several of the identified genes correspond to known drug targets, while others represent potential candidates for further investigation. These results illustrate how integrating molecular alterations with network-based modeling can help stratify glioblastoma tumors and uncover molecular mechanisms that may guide the development of more personalized therapeutic strategies.

The challenge of stratifying patients for combination therapy is both technically demanding and clinically crucial. In previous work, we introduced a multi-objective optimization framework for identifying optimally synergistic combination protocols that are robust to competing definitions of additivity. This manuscript extends this methodology to quantify how inter-individual variability in drug sensitivity influences the combination doses that optimally balance the competing objectives of synergy of efficacy and synergy of potency (a proxy measure of toxicity). For this methodology, we introduce a voxel-based stratification approach to characterize individuals (model parameterizations) into subgroups based on sensitivity to each drug as a monotherapy and in combination. As a case study, we apply the method to a preclinical dataset of murine response to the combination of an immune checkpoint inhibitor and an antiangiogenic agent. We demonstrate that the algorithm can quantify how the robustly optimal combination therapies vary across different treatment response subgroups and how the algorithm can identify subpopulations for which no meaningfully efficacious combination exists. As applying the methodology requires knowledge of specific parameter values for which measurable biomarkers may be unavailable, we also propose an initiation protocol that permits identification of the parameters necessary to place an individual in a subgroup. This methodology is a step in the direction of determining the right combination therapy for a subgroup and finding the right subgroup for an existing therapy.

Metabolism plays an essential role in cellular processes: development, growth, differentiation, and determination of cell identity. Understanding how metabolic processes dynamically change across cell types, stages, and environmental conditions is crucial for studying developmental biology, aging, and disease progression. Genome-wide metabolic models (GEMs) are a powerful network-based tool for studying these processes by integrating omics data to model context-specific metabolism. However, current approaches, such as Flux Balance Analysis (FBA), have limitations in addressing the dynamic nature of metabolism across developmental stages at a sample-specific resolution. To address this, we introduce a novel network-based method for analyzing cell and stage specific metabolic flow using directed and weighted metabolic networks that account for sample-specific transcriptomic data. We apply this method to study ovarian follicle development, providing a deeper understanding of intra-cellular metabolic processes, identifying key metabolites, enzymes, and potential markers for follicular maturation, important for IVF. By incorporating biologically meaningful data, this approach bridges the gap between theoretical metabolic network models (GEMs) and experimental observations, offering a systems-level view of metabolic dynamics in developmental and understudied contexts.

The aryl hydrocarbon receptor (AhR) is a ligand-activated transcription factor involved in xenobiotic sensing, as well as development, immunity, and tissue homeostasis. AhR signaling can proceed through a canonical and non-canonical pathway; the present study focuses on the canonical pathway. While ligand-dependent differences in binding affinities and direct ligand degradation kinetics are well known, and subtle differences in ligand binding can shape downstream signaling, it is still unclear which biochemical reaction steps within the canonical pathway are responsible for distinct ligand-specific transcriptional responses. Here, we developed a mechanistic ordinary differential equation model of the canonical AhR pathway. We calibrated the model to time-resolved qPCR measurements of Cyp1a1 and Ahrr mRNA in mouse bone-marrow-derived macrophages exposed to structurally diverse, environmentally relevant ligands with known immunomodulatory activity (3-methylcholanthrene, indolo[3,2-b]carbazole, and bisphenol A) using global optimization under a heteroskedastic likelihood. To dissect ligand specificity, we evaluated 528 candidate models that allow one or two ligand-involving reaction rate constants to vary. Akaike-based model selection reveals a dominant dynamical regime governed by promoter occupancy and target-gene mRNA synthesis, indicating that ligand-specific transcriptional responses are primarily encoded at the level of transcriptional regulation rather than upstream signaling events. The resulting model is made available in SBML and PEtab formats for reproducibility, and to enable further research into whether ligand-specific effects are conserved or rewired across cell types.

Forecasting how microbiome-host ecosystems evolve through time simultaneously at the compositional and functional level remains a central challenge in biology. While dynamical systems models (DSMs) can infer and predict community composition from longitudinal abundance data, and constraint-based metabolic models (CBMMs) can estimate metabolic fluxes from genome-scale reconstructions, no existing framework unifies these approaches to generate mechanistically grounded, time-resolved forecasts of both microbial abundances and metabolite dynamics from ecological data alone. Here, we introduce the Dynamical Systems Constrained Metabolic Modeling (DySCoMeMo) framework, a new hybrid computational pipeline that integrates ecological DSMs with CBMMs to predict temporal dynamics of biomass and metabolites across microbial communities and hosts. DySCoMeMo leverages parameters inferred from application of DSMs to microbiome time series data to constrain metabolic modeling over time, thereby bridging ecological interaction networks with genome-scale metabolic modeling. DySCoMeMo is able to predict future community and metabolite dynamics in vitro with accuracy that is superior or on-par compared to that achieved with established methods that require actual microbial abundances and/or metabolites data for metabolite network inference or for estimating the per-microbe contribution to the extracellular metabolic pool. DySCoMeMo also generalizes to in vivo data as it is capable of accurately forecasting microbial and metabolite dynamics in response to dietary perturbations even when host metabolism is included. Finally, DySCoMeMo uniquely enables the identification of keystone species by quantifying their contributions to sustaining metabolic environments. Together, our work establishes a generalizable, mechanistically grounded framework for time-resolved forecasting of microbiome-host microbial and metabolic dynamics, bridging ecological interaction inference with genome-scale metabolism of communities.

Migration of leukocytes in the context of immune homeostasis or inflammatory diseases is regulated by activation of chemokine receptors by chemokine ligands. To elucidate how these interactions give rise to cell migration, we mapped the chemokine-stimulated signal transduction network in monocytic THP-1 cells. Global phosphoproteomics revealed 630 time-resolved changes in phosphorylated proteins downstream of the chemokine receptor CCR2. We used the "PHONEMeS" network modeling algorithm to generate the most parsimonious signal transduction network consistent with the observed protein phosphorylation data. The CCR2 signaling network is highly divergent, acting via multiple branches to regulate proteins required for cell migration. We validated this model using kinase inhibitors targeting different branches of the network and successfully blocked chemokine-stimulated cell migration. Thus, chemotaxis is an emergent property resulting from an integrated cellular response to divergent signaling pathways. This paradigm suggests that physiological regulation or pharmacological blockade of chemokine-driven inflammation could potentially be achieved by inhibiting any of the divergent pathways within the network.

Time-resolved measurements are central to calibrating mechanistic dynamical models, but current inference frameworks typically assume that reported measurement times are exact. In practice, actual sampling times may deviate from reported times because of sample-handling delays, imper-fect synchronization, or reporting errors. Here, we present a Bayesian framework for parameter inference in ordinary differential equation models that explicitly accounts for uncertainty in measurement times. We formulate latent measurement times as random variables and derive a joint and marginalized posterior. To compute the marginal likelihood efficiently, we augment the original dynamical system with additional state variables that evaluate the required integrals during numerical simulation. This reduces the dimensionality of the estimation problems and allows for efficient and reliable Markov chain Monte Carlo sampling. Across synthetic examples and a published model of carotenoid cleavage in Arabidopsis thaliana, neglecting time uncertainty led to biased estimates and overconfident uncertainty quantification, whereas the proposed marginalized formulation recovered reliable parameter estimates while substantially improving sampling efficiency and scalability. These results identify measurement time uncertainty as an important source of variability in dynamic modeling and establish posterior marginalization as a practical strategy for robust mechanistic inference.

The long-held view that acetate, one of the main fermentation by-products of Escherichia coli, is toxic to microbial growth is currently challenged. Here, we demonstrate that acetate promotes E. coli adaptation to nutrient changes by accelerating growth resumption, with as little as 250 {micro}M acetate being sufficient to shorten the lag phase by several hours. Acetate was found to be consumed via acetyl-CoA synthetase very early after the nutrient change. Transcriptomics, metabolomics and 13C-isotope labeling experiments show that acetate replenishes metabolic pools in the tricarboxylic acid cycle and upper glycolysis. Single-cell analyses reveal that acetate increases the adaptation speed of individual cells switching to the new nutrient. We conclude that the reuse of excreted acetate by E. coli facilitates metabolic adaptation by transiently replenishing central metabolite pools. This work identifies an unexpected role of acetate in the nutritional adaptation of E. coli, providing new insights into the physiological relevance of overflow metabolism. HighlightsO_LIAcetate facilitates E. coli adaptation from one nutrient to another. C_LIO_LILess than 250 {micro}M acetate is sufficient to halve lag times. C_LIO_LIAcetate helps replenish metabolite pools in central carbon metabolism. C_LIO_LIAcetate excretion is an adaptative strategy to overcome resource fluctuations. C_LI

Reconstructing dynamic immune cell state transitions from single-cell transcriptomic data requires coordinated analytical strategies that capture both phenotypic progression and underlying regulatory programs. This protocol describes a step-by-step computational workflow for analyzing human tumor-infiltrating T cells using the sequential application of dimensionality reduction, pseudotime trajectory inference, regulon activity analysis, and transcription factor-transcription factor network reconstruction. The workflow outlines data preprocessing and quality control, trajectory rooting and parameter selection, branch-specific differential analysis, and the integration of regulon inference to contextualize transcriptional programs along inferred trajectories. Regulon-based TF-TF network reconstruction is used as a downstream interpretive layer to identify regulatory modules associated with distinct cell-state transitions. Publicly available at GitHub repository https://github.com/rogercasalsfr/immuno-trajectory-grn-integrative-workflow, this protocol emphasizes practical considerations including parameter sensitivity, trajectory robustness, and consistency between phenotypic and regulatory outputs. The protocol supports reproducible analysis and interpretation of immune cell dynamics in human tumor microenvironment studies using single-cell RNA sequencing data.